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ABSTRACT 

We have performed nonlinear ideal magnetohydrodynamic simulations of the 
long term evolution of a magnetized low-density "bubble" plasma formed by 
a radio galaxy in a stratified cluster medium. It is found that about 3.5% of 
the initial magnetic energy remains in the bubble after ~ 8 x 10 9 years, and 
the initial magnetic bubble expansion is adiabatic. The bubble can survive for 
at least 8 x 10 9 years due to the stabilizing effect of the bubble magnetic field 
on Rayleigh- Taylor and Kelvin-Holmholtz instabilities, possibly accounting for 
"ghost cavities" as observed in Perseus-A. A filament structure spanning about 
500 kpc is formed along the path of bubble motion. The mean value of the 
magnetic field inside this structure is ~ 0.88 fiG at ~ 8 x 10 9 years. Finally, the 
initial bubble momentum and rotation have limited influence on the long term 
evolution of the bubble. 

Subject headings: galaxies: jets — magnetic fields — MHD — methods: numerical 



1. Introduction 



An unsolved problem in active galactic nuclei (AGN) feedback on clusters is how to 
account for the the morphology and stabil ity of buoyant bubbles and th eir interactions with 
the ambient intracluster medium (ICM) (IMcNamara fc Nulsenl 120071 ). iFabian et al.l (120061 ) 
showed that in the Perseus cluster, such bubbles can stay intact far from cluster centers where 
they were inflated by AGN jets. However, studies of kinetic-energy dominated jets in the 
purely hydrodynamic limit did not explain the observed long term persistence of the buoyant 
bubbles, which are prone to Rayleigh- Taylor (RTI) and Kelvin-Holmhotz (K HI) instabilities 
and fr a gment entirely within 10 Myr, contrar y to observations [but see iReynolds et al. 
f l2005h : Ipizzolato fc Sokerl J2006h : lGardii3 (b()07t )]. 
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Appreciable magnetic energy has be e n observed in both c luster and radio lobe plasmas 
f lOwen et al.l I2OO0I : iKronberg etaD l200ll : ICroston et all |2005| ). The magnetic fields could 
play a vital role in the dynamics of the rising bubble, as shown by a series of 2D magne- 



tohydro dynamic (MHD) studies for bubbles o f several 10 8 years (IBriiggen &: Kaiser! 12001 



Robinson et al. 2004 ; Jones fc De Young 2005 ). Stone fc Gardiner ( 20081 ) showed that uni- 



form strong magnetic fields do not suppress RTI completely, but sheared fields at the bubble 
interface quench the instability. This impli es that not only the field strength but also the 
configuration matters for bubble stability. iRuszkowski et al.l (120071 ) did a comprehensive 
study of the influence of different magnetic field configurations upon the stability of a rising 
bubble. They found that the internal bubble helical magnetic field moderately stabilizes 
the bubble, however the bubble has an initial plasma parame ter (3 = (2nT/B 2 ) ^> 1, still 
in thermal energy-dominated regime. Recently, iLi et al.l (120061 ) proposed a differen t bubble 



magne tic field configuration for jet/lobes. Here we adopt the field configuration of ILi et al. 



(120061 ) and focus on the late stage of the evolution of magnetically dominated bubbles in 3D. 
We show that this spheromak-like magnetic field configuration strongly stabilizes the insta- 
bilities and prevents the bubble from breaking up, possibly forming the intact but detached 
"ghost cavity" observed in systems such as Perseus-A. This Letter is organized as follows. 
In Sec. El we describe the problem setup. The simulation results and discussions are given 
in Sec. |3j 



2. Problem Set Up 



The background ICM is assumed to be hydrostatic and isothermal. The density and 
pressure profile is given as p = p = [1 + (R/R C ) 2 ]~ K , where the parameters k and R c are 
taken to be 1.0 and 4.0 respectively. The fi xed yet distributed gra vitational field — V^(-R) 
dominated by the dark matter is assumed (INakamura et al.ll2006l ). We simulate only the 
phase after the AGN has inflated a low-density bubble which is in pressure equilibrium 
with the ICM initially. The static magnetized bubble initially lies at z 6 ) = (0,0,5) 

with density p b = 0.1, spherical radius r d = 2 for the density profile and r b = 1 for the 
magnetic configuration. The specific heat 7 inside and outside the bubble is taken to be 5/3. 
Physical quantities are normalized by the characteristic system length scale Ro = 25 kpc, 
density po = 1-67 x 10~ 26 g cm -3 , and velocity Vq = 6.2 x 10 7 cm s _1 . The initial sound 
speed, C s | i=0 = 7 1 / 2 ~ 1.29, is constant throughout the computational domain. Other 
quantities are normalized as: time t = 1 gives Rq/Vq = 38.6 Myr, magnetic field B = 1 
1.13 x lO" 5 G and energ y E = 1 gives p n VnRl = 2.71 x 10 58 ergs. 



2U/2 



gives (2p Vo 

The magnetic field setup of the bubble follows ILi et al.l (|2006l ) with the ratio of toroidal to 
poloidal bubble field a taken to be a/10, which corresponds to a minimum initial Lorentz 
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force. The latter is reasonable since we want to mimic the situation that the magnetic 
bubble has substantially relaxed and detached from the jet tip. The initial total magnetic 
energy is around 2.0 x 10 59 ergs. The computational domain is taken to be \x\ < 16, 
\y\ < 16, and < z < 32 corresponding to a (800 kpc) 3 box in the actual length scales. The 
numerical resolution used here is 400 3 , where the grid points are assigned uniformly in every 
direction. A cell 5x corresponds to 2.0 kpc. We use "outflow" boundary conditions at every 
boundary except the boundary at z = 0, where we use "reflecting" boundary conditions, 
which guarantee that the energy flux through this boundary is zero. 



3. Results & Discussion 

Energy and density evolution reveal the different stages of the bubble. Figure W^top) 
presents the time evolution of various energies at the early stage, in which the gravita- 
tional energy E g = J pipdV and the internal energy E? = fp/{l — ^-)dV, where dV is 
the infinitesimal volume and the integral is over the entire computation domain, have been 
reduced by their initial values E gj0 and Et,o, respectively. The total energy is defined as 
-Etotai = E m + Ek + Et + E 9 , where the kinetic energy is defined as Ek = J l/2pv 2 dV and the 
magnetic energy is defined as E m = J B 2 /2dV. At the early stage (t < 5), the kinetic energy 
increases since the bubble accelerates upward due to buoyancy, and the magnetic energy de- 
creases due to work don e in expansion from th e weak initial Lorentz force and conversion to 
shock and wave energy jNakamura et alJbood ). The gravitational energy increases because 



there is a net outward mass flow in the axial direction. The passage of the shock wave heats 
and compresses the ICM and alters its pressure gradient. The total energy £7 to ta,i is almost 
constant before t ~ 10, at which time the shock and wavefront reach the boundary. Conser- 
vation of the total energy for t < 10 is not strictly satisfied due to numerical diffusion. The 
initial bubble expansion (t < 20) is approximately adiabatic, which gives E = Vp oc V 1 ^"', 
where V is the bubble volume. Therefore EijE\ = (V1/V2) 2 / 3 , where the subscript 1 and 2 
indicate the initial and final state, respectively. At t — 20, the radius of the bubble has in- 
creased from 2 to ~ 5, which gives E 2 /Ei ~ 16%, roughly matching the amount of magnetic 
energy in the bubble [20% from Fig. [T](iop)]. For later times [Fig. [T]( bottom)], fitting with 
E m (t) / E m fl = exp(— t/rdi S ), we have: (1) a fast dissipation stage (t < 10) when magnetic 
energy dissipation time Tdi S ~ 11; (2) a slow dissipation stage (t > 10) when Tdi S ~ 114. 
Both times are much smaller than the numerical dissipation time r res at the corresponding 
stages (see discussions at the end of this section). The kinetic energy is oscillating while it is 
decaying slowly, which can be understood as follows. Gravity pulls the bubble down to the 
denser ICM, compressing the bubble and causing an increase in magnetic energy (since the 
magnetic flux is nearly conserved). This results in the magnetic energy oscillating in phase 
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with the kinetic energy with period ~ 90 (3.5 x 10 9 yr). At t — 200 after several periods of 
decaying oscillations, about 3.5% of the initial magnetic energy remains. 

The magnetic field suppresses instabilities and therefore the bubble remains intact 
longer. Figure [2] presents the typical density distribution (logarithmic scale) in 2D x—z 
slices at y — at different times (t = 0, 7.5, 20, 50, 100, 125). The white solid contour lines 
indicate contours of constant magnetic field strength \B\. Consistent with Fig. ^bottom), 
after t > 50, the bubble undergoes a slowly decaying oscillation between z ~ 19 and z ~ 24. 
We have also performed simulations of an unmagnetized bubble and found that it disinte- 
grates after t = 20, whereas the magnetic bubble still clearly differentiates itself from the 
ambient medium at t — 200 (Fig. [2]). The formation of an "umbrella" or a thin protec- 



tive m agnetic layer on the bubble working surface suppresses instabilities (jRuszkowski et al. 



20071 ). The stronger magnetic field is also found to move the bubble faster and push the 
bubble farther away from its initial position. The position of the bubble top as a function of 
time is shown in Fig. [31 This position is calculated by plotting the axial profile of the density 
along the line of (x, y) = (0, 0) and finding the location of the first density jump from the top 
of the domain. Comparing the unmagnetized run (dash line) with the magnetized run (both 
without initial momentum/rotation), we can see that the rising speed increases from 0.36 to 
0.6. Note that the rough estimate of the terminal speed v t based on the initial gravitational 
acceleration g and size rj, both position-dependent, is v t ~ 4/3v / 2qr v; ~ 2, assuming force 



balan ce between buoyancy and viscous drag force at the final stage (IMcNamara fc Nulsen 



2003). 



We investigate the effe ct of the magnetic field on bubble stability in more detail. Follow- 



ing [N^kamui^eLalJ (120071 ). we first study KHI of the bubble. At both sides of the bubble, the 
magnetic field is almost uniform, not twisted like at the top of the bubble. From linear anal- 
ysis, the instability criterion for nonaxisymmetric KHI surface modes is (jHardee fc Rosen 



20021 ) : AV > V As = [(p b + p e )/(Anp b p e )(B 2 + B 2 )] 1/2 , where AV = \V b - V e | is the veloc- 
ity shear and Vks is the surface Alfven speed. The subscripts b and e indicate the bubble 
an d external medium, res pectively. The nonaxisymmetric body modes would be important 



if (jHardee fc Rosenlll999l ): (1) V b > Vf, in which V/ is the fast magnetosonic speed; or (2) 
C s Va/ (Cg+V^) 1 ^ 2 < V b < V s , in which Va is the Alfven speed and V s is the slow magnetosonic 
speed. The definitions of V u are: V u = {1/2{C S 2 + Vf ± [(C 2 S + Vj) 2 ~4:C 2 S V 2 cos 2 9} 1 / 2 }} 1 ' 2 . 
Since we focus on the x-direction only, Vf cos 2 9 are taken to be B 2 / p. FigureH] (top) displays 
the transverse distribution of the bulk flow speed V = (v 2 + v 2 + v 2 ) 1 ^ 2 , the density p, and 
the magnetic field strength B at t — 7.5 at z = 8.0. A distinct velocity shear is identified 
at x ~ 3.5. Across this shear, AV ~ 0.75 and Va s ~ 0.85. The inequality AV < Va s holds. 
This means that the KHI surface modes are completely suppressed. FigureH] (bottom) dis- 
plays the transverse distribution of V, V/, and V s at t = 7.5 at z — 8.0 . We see that the 
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bulk flow V lies between V) and V^. The inequality V s < V < Vf is satisfied in the body of 
the buoyant bubble. This rules out the KH1 body modes as well. 

RTI could also be suppressed by the magnetic field. For the idealized case of two 
conducting fluids separated by a contact discontinuity with a uniform magne tic field parallel 
to the interface undergoing constant acceleration g, I Chandr asekharl (119611 ) demonstrated 
that RTI on a scale L parallel to the field requires B < B c = [Lg(ph — pi)] 1 ^ 2 , in which ph 
and pi are the densities in the heavy and light fluids respectively. Modes perpendicular to 
the field are unaffected. At the top of the bubble (z ~ 11.6) (Fig. [5]), ph and pi are found to 
be 0.15 and 0.06, respectively, and the gravitational acceleration g is calculated to be 0.95. 
If L is chosen to be the computation domain size, the maximum possible mode wavelength 
in the simulations, then the critical magnetic field strength is B c ~ 1.6, which is larger 
than the magnetic field (~ 0.8) at that location. This means that only part of the parallel 
modes are suppressed regardless of the perpendicul ar modes. This is not c onsistent with 
the simulation results. However, as pointed out in IStone fc Gardinerl (120081 ). the twisting 
nature of the bubble field at the top of the bubble introduces a current sheet at the surface, 
and the changes in the direction of the field at the interface must be on very small scales to 
inhibit the interchange modes. A more detailed study of this effect is beyond the scope of 
this paper and will be the subject of future study. 

Our simulation also provides one possible explanation of the morphology and origin 
of the large s cale magnetic fi e ld an d the generation of "ghost cavities" observed in many 
clusters. Like iRobinson et al.l (120041 ) . a magnetized high-density tail remains as the bubble 
rises (Fig. [2]). The tails have an elongated morphology (Fig. [2]), resembling Ha filaments, 
which is found to indicate the history of the rising bubble. Interest ingly, the magnet i c field 
also helps stabilize this "filament" as it is still visible at t = 200. iNipoti fc Binneyl (120041 ) 
has argued that thermal conduction has to be strongl y suppressed in the ICM ; otherwise, 
such cold filaments would be rapidly evaporated. As in iRuszkowski et al.l (120071 ). our results 
provide a possibility that thermal conduction may be locally weaker in the bubble wake, thus 
preventing or slowing down filament evaporation. The magnetic field is distributed between 
z ~ 2.5 and z ~ 25 at t = 200 with peak value around 0.156 (~ 1.8 uG ) or mean value 
around 0.078 (~ 0.88 u G), close to the estimates of wider cluster fields (jCarrili fc Taylor 



20021 ; iTaylor et al.ll2002l ). This large-scale (> 500 kpc) magnetic field structure also mimics 
the morphology of the second class ("Phoenix ") of "radio relics" elongated from the cluster 
center to the periphery observed in A 115 (jGovoni et al.l l200ll ). The simulation results 
mean that the magnetic bubble rising from ~ 125 kpc with an initial magnetic energy of 
2 x 10 59 ergs, will spread magnetic fields between 62.5 kpc < z < 600 kpc and keep a 
magnetized bubble between 475 kpc and 600 kpc for at least ~ 8 x 10 9 yr. Therefore one 
could reproduce the intact but detached "ghost cavity" observed in some systems such as 
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Perseus-A. 

Real radio lobes possess complex, jet-driven flows, different from the static bubble 
that we use for our initial conditions. We study the effects of additional variations in the 
initial bubble by having: (1) an initial uniform injection speed v- m j or (2) an initial uni- 



form rotation uj. iNakamura et al.l (120081 ) reported that the radio lobe gains a momentum 
of finj ~ 0.91C s0 = 1.174, and the rotation speed u at the edge of the bubble reaches 
~ l-hCs^/r^ = 0.9675 after jet-driven inflation and interaction with the ambient ICM, al- 
though neither are not uniform inside the bubble. Here we take these two values as our 
initial values for t^j and uj and idealize both to be uniform throughout the bubble. Fig. [3] 
presents the time evolution of the position of the bubble top with initial uniform injection 
velocity t> in j = 1.174 (long dash) and initial uniform rotation speed uj = 0.9675 (dash dot). 
It is seen that these two variables only have some influence on the evolution of the bubble 
during the early stage. The bubble with an initial rotation has a reduced buoyancy resulting 
from rotation-driven instabilities that in turn lead to a reduced pressure differential between 
the bubble and the ICM. The evolution, and structure of the bubble in the later times are, 
however, essentially the same. 

The effects of numerical diffusion are estimated as follows. If we fit the net toroidal 
magnetic flux ip t = J B y dS (only positive B y is selected) as ^t{t) /ipt(t — 0) = exp(— t/r Tes ), 
the "resistive" dissipation time due to numerical diffusion is r ros ~ 38 for t < 10, r res = 138 
for 10 < t < 50, and r res = 507 for 50 < t < 100. Therefore the numerical diffusion is 
not important on the time scales of interest (t < 200). Increasing numerical resolution can 
further reduce the numerical diffusion. 



The oldest bubbles known so far are a few hundred Myr old (IMcNamara fc Nulsen 



20071 ). which are much younger than the long-sustained (~ a Hubble time) magnetized 
bubbles reported in this Letter. The simulations show that at later times t > 50, the density 
and X-ray luminosity contrast between the bubble and ICM plasma at the periphery of the 
cluster is small, which makes the direct X-ray observation of the bubble faraway from the 
parent galaxy very difficult. But the radio bubble at the borders of the clusters may still 
be detectable through: (1) Faraday rotation measurement of the outskirts of the cluster if 
an outside synchrotron source exists; (2) or radio observation from re-energization of fossil 
radio bubble plasma by shocks or mergers. It is highly possible for the bubbles to pile up in 
the outer atmospheres of the clusters in deep Chandra observation, if those observations are 
technically viable. 

It is conceivable that part of the initial magnetic energy of the bubble accelerates cosmic 
rays. Throughout the lifetime of a cluster, there could be multiple AGNs injecting jets/lobes 
into the ICM. Both the cosmic rays and the remaining magnetic fields (distributed over 
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large scales) could provide the energy sources for pheno mena such as radio relics and radio 



haloes that have been observed for a number of clusters (IFerrari et al.l 120081 ) 
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Fig. 1. — Time evolution of various energies. Top: short term various energies evolution 
t < 20, thin solid: magnetic Energy E m , thin dash: kinetic energy Ek, thin dash dot: internal 
energy E?, thick dash: gravitational energy E g and thick solid: total energy E to \. Bottom: 
long term time evolution of magnetic energy (in logarithmic scale) E m (solid) and kinetic 
energy Ek (dash), t — 1 represents 38.6 Myr. 
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Fig. 2.— (color) Density (in logarithmic scale) in the x-z plane as a function of time 
(a = y/lO). The white solid contour lines indicate the magnetic field strength \B\. t = 0, 
\B\ G [0,1.976]; t = 7.5, \B\ G [0,0.830]; t = 20, \B\ e [0,0.400]; t = 50, \B\ E [0,0.328]; 
t = 100, \B\ E [0,0.159]; t = 125, \B\ G [0,0.156]. \B\ = 1 represents 11.3 fiG. The number 
of the contour levels are all 5. 
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Fig. 3. — Position of the bubble top v.s. time for four different runs. The solid, long dash 
and dash dot lines represent magnetized bubble runs without initial momentum rotation, 
with f in j = 1.174 and with uo = 0.9675 respectively. The short dash line denotes the purely 
hydrodynamic run. 
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Fig. 4. — Transverse profiles in the x-direction of several quantities at t — 7.5 on z = 8.0. 
The bulk speed V, the density p, and the magnetic field strength B are shown for inspecting 
the KHI surface modes (top). The bulk speed V, the superfast magnetosonic speed V), and 
the slow magnetosonic speed V s are shown for inspecting the KHI body modes (bottom). 
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Fig. 5. — Axial profiles in the ^-direction of several quantities at t — 7.5 with (x, y) = (0, 0). 
The density p and the magnetic field strength B are shown for inspecting the RTI modes. 



